package ua.lviv.franko.solvers.stressanalysis.twoDimention;

import java.util.ArrayList;
import java.util.Iterator;
import java.util.List;

import ua.lviv.franko.integration.Gauss2D;
import ua.lviv.franko.matrix.RowType;
import ua.lviv.franko.matrix.hash.lu.HashMatrixLUStyle;
import ua.lviv.franko.mlsFunctions.twoDimension.BaseFunction;
import ua.lviv.franko.mlsFunctions.twoDimension.Particle;
import ua.lviv.franko.mlsFunctions.twoDimension.Values;
import ua.lviv.franko.mlsFunctions.twoDimension.weightFunctions.IWeightFunction2D;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.fem.AbstractStressElement;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.savestructures.EquationParameters;

public class CoupledSolver {
	public static final double					EPS						= 0.000000000001;
	private static final int					NUMBER_OF_UNKNOWNS		= 2;
	private static final int					QADRATURE_POINTS_NUMBER	= 19;
	// general data
	protected Gauss2D							gau						= new Gauss2D(QADRATURE_POINTS_NUMBER);
	protected EquationParameters				parameters;
	protected HashMatrixLUStyle						mainMatrix;
	protected double[]							b;

	// data for base functions
	protected ArrayList<Particle>				particles;
	protected ArrayList<BaseFunction>			baseFunctions;
	protected ArrayList<AbstractStressElement>	elementsList;
	int											N;

	public CoupledSolver(EquationParameters parameters, ArrayList<Particle> particles, ArrayList<AbstractStressElement> elementsList, int n, int PolinomialDegree,
			IWeightFunction2D weight) {
		super();
		this.parameters = parameters;
		this.particles = particles;
		this.elementsList = elementsList;
		N = n;

		createBaseFunctions(PolinomialDegree, weight);
	}

	public void generate() {
		mainMatrix = new HashMatrixLUStyle(RowType.ARRAY_TYPE, (N + particles.size()) * NUMBER_OF_UNKNOWNS);
		b = new double[(N + particles.size()) * NUMBER_OF_UNKNOWNS];

		// first build matrix for finite elements
		System.err.println("start fem matrix building");
		for (Iterator<AbstractStressElement> iterator = elementsList.iterator(); iterator.hasNext();) {
			AbstractStressElement element = iterator.next();

			ElementMatrixForStress data = element.generateAll();
			List<Integer> pos = data.getNodeList();
			for (int i = 0; i < pos.size(); i++) {
				for (int j = 0; j < pos.size(); j++) {
					mainMatrix.setValue(pos.get(i) * 2, pos.get(j) * 2, mainMatrix.getValue(pos.get(i) * 2, pos.get(j) * 2) + data.matrix[2 * i][2 * j]);

					mainMatrix.setValue(pos.get(i) * 2, pos.get(j) * 2 + 1, mainMatrix.getValue(pos.get(i) * 2, pos.get(j) * 2 + 1) + data.matrix[2 * i][2 * j + 1]);
					mainMatrix.setValue(pos.get(i) * 2 + 1, pos.get(j) * 2, mainMatrix.getValue(pos.get(i) * 2 + 1, pos.get(j) * 2) + data.matrix[2 * i + 1][2 * j]);
					mainMatrix.setValue(pos.get(i) * 2 + 1, pos.get(j) * 2 + 1, mainMatrix.getValue(pos.get(i) * 2 + 1, pos.get(j) * 2 + 1) + data.matrix[2 * i + 1][2 * j + 1]);
				}
				b[pos.get(i) * 2] += data.vector[i * 2];
				b[pos.get(i) * 2 + 1] += data.vector[i * 2 + 1];
			}
		}
		System.err.println("finish fem matrix building");
		// now build matrix for efg
		System.err.println("start efg matrix building");
		int femDisplacement = N * NUMBER_OF_UNKNOWNS;
		for (int i = 0; i < b.length / NUMBER_OF_UNKNOWNS; i++) {
			BaseFunction baseFunction = baseFunctions.get(i);
			double[][] vals = this.gaussQuadratureIntegration(baseFunction.getXi().getXi() - baseFunction.getXi().getRho(), baseFunction.getXi().getXi()
					+ baseFunction.getXi().getRho(), baseFunction.getXi().getYi() - baseFunction.getXi().getRho(), baseFunction.getXi().getYi() + baseFunction.getXi().getRho(),
					baseFunction, femDisplacement);

			int pos1 = particles.indexOf(baseFunction.getXi());
			for (int j = 0; j < vals.length; j++) {
				int pos2 = particles.indexOf(baseFunction.getNodeList().get(j));
				mainMatrix.setValue(pos1 * NUMBER_OF_UNKNOWNS + femDisplacement, pos2 * NUMBER_OF_UNKNOWNS + femDisplacement, vals[j][0]);
				mainMatrix.setValue(pos1 * NUMBER_OF_UNKNOWNS + femDisplacement, pos2 * NUMBER_OF_UNKNOWNS + 1 + femDisplacement, vals[j][1]);
				mainMatrix.setValue(pos1 * NUMBER_OF_UNKNOWNS + 1 + femDisplacement, pos2 * NUMBER_OF_UNKNOWNS + femDisplacement, vals[j][2]);
				mainMatrix.setValue(pos1 * NUMBER_OF_UNKNOWNS + 1 + femDisplacement, pos2 * NUMBER_OF_UNKNOWNS + 1 + femDisplacement, vals[j][3]);
			}
		}
		System.err.println("finish efg matrix building");
		// print global matrix
		System.out.println("-----------------");
		for (int i = 0; i < b.length; i++) {
			for (int j = 0; j < b.length; j++) {
				System.out.print(mainMatrix.getValue(i, j) + "\t");
			}
			System.out.print("\n");
			// System.out.println(b[i]);
		}
		System.out.println("-----------------");

	}

	private void createBaseFunctions(int PolinomialDegree, IWeightFunction2D weight) {
		this.baseFunctions = new ArrayList<BaseFunction>();
		for (Iterator<Particle> iterator = particles.iterator(); iterator.hasNext();) {
			Particle particle = iterator.next();
			BaseFunction bf = new BaseFunction(PolinomialDegree, particle, weight);
			bf.setNodeList(this.getParticleSupport(particle, weight));
			this.baseFunctions.add(bf);
			// bf.display();
		}
	}

	private double[][] gaussQuadratureIntegration(double ax, double bx, double ay, double by, BaseFunction baseFunction, int femDisplacement) {
		// first column - baseFunction*v; second column - v*baseFunction
		double[][] result = new double[baseFunction.getNodeList().size()][NUMBER_OF_UNKNOWNS * NUMBER_OF_UNKNOWNS];
		// its a position of main node - node of base function
		int pos = baseFunction.getNodeList().indexOf(baseFunction.getXi());

		double f1 = 0.0;
		double f2 = 0.0;

		for (int j = 0; j < QADRATURE_POINTS_NUMBER; j++) {
			double temp1 = 0.0;
			double temp2 = 0.0;
			for (int i = 0; i < QADRATURE_POINTS_NUMBER; i++) {
				double x = (ax + bx) / 2 + (bx - ax) / 2 * gau.getX(i);
				double y = (ay + by) / 2 + (by - ay) / 2 * gau.getX(j);
				Values[] values = baseFunction.calculateAllFunctions(x, y);
				double k1 = parameters.getK1();
				double k2 = parameters.getK2();
				double k3 = parameters.getK3();

				for (int k = 0; k < values.length; k++) {
					result[k][0] += gau.getW(j) * gau.getW(i) * (values[pos].valX * values[k].valX * k1 + values[pos].valY * values[k].valY * k3);
					result[k][1] += gau.getW(j) * gau.getW(i) * (values[pos].valY * values[k].valX * k2 + values[pos].valX * values[k].valY * k3);
					result[k][2] += gau.getW(j) * gau.getW(i) * (values[pos].valY * values[k].valX * k2 + values[pos].valX * values[k].valY * k3);
					result[k][3] += gau.getW(j) * gau.getW(i) * (values[pos].valX * values[k].valX * k3 + values[pos].valY * values[k].valY * k1);

				}

				temp1 += gau.getW(i) * values[pos].val * parameters.getF1().calculate(x, y);
				temp2 += gau.getW(i) * values[pos].val * parameters.getF1().calculate(x, y);
				// temp += gau.getW(i) * coefFunc.calculate(x, y);
			}
			f1 += temp1 * gau.getW(j);
			f2 += temp2 * gau.getW(j);
		}
		f1 *= (bx - ax) * (by - ay) / 4.0;
		f2 *= (bx - ax) * (by - ay) / 4.0;
		for (int i = 0; i < result.length; i++) {
			result[i][0] *= (bx - ax) * (by - ay) / 4.0;
			result[i][1] *= (bx - ax) * (by - ay) / 4.0;
		}
		// result[pos][1] = f;
		b[particles.indexOf(baseFunction.getXi()) * NUMBER_OF_UNKNOWNS + femDisplacement] = f1;
		b[particles.indexOf(baseFunction.getXi()) * NUMBER_OF_UNKNOWNS + 1 + femDisplacement] = f2;
		return result;
	}

	private ArrayList<Particle> getParticleSupport(Particle p, IWeightFunction2D weight) {
		ArrayList<Particle> list = new ArrayList<Particle>();
		for (Iterator<Particle> iterator = particles.iterator(); iterator.hasNext();) {
			Particle particle = (Particle) iterator.next();
			if (weight.intersectedSupport(p, particle)) {
				list.add(particle);
			}
		}
		return list;
	}

}
